#include "HY1C_out.h"
/* ========================================================================
 *  Procedures create, read, and write HDF mscalibration pixel data
 *
 * int crosscal_append(char *xcalfile, mscalstr calstr) - creates HDF, if
 *			it does not exist already, and appends data
 * int crosscal_read(char *xcalfile, int32_t subsmpl, mscalstr *calstr) - 
 *			reads subsampled pixel cross-calibration data and
 * 			allocates memory and data to the calstr structure
 * int crosscal_add(char *xcalfile, int32_t subsmpl, mscalstr *calstr, int32_t *npixs, int32_t *ngranuls) 
 *                      reads cross-calibration data from an hdf file,
 *          		subsamples if needed, and adds pixel and granule
 *          		data to an already existing calstr structure
 *
 *     Programmer     Organization      Date       Description of change
 *   --------------   ------------    --------     ---------------------
 *   Ewa Kwiatkowska  SAIC           10 September 2003    Original development
 *   Joel Gales       Futuretech     24 October   2012    Add ds_id.fftype
 *
 * ======================================================================== */


#include <stdio.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <stdlib.h>
#include <fcntl.h>
#include <string.h>
#include <errno.h>
#include <unistd.h>
#include <libgen.h>
#include <math.h>

#include "passthebuck.h"
#include "l2prod_struc.h"
#include "filehandle.h"
#include "hdf_utils.h"
#include "mfhdf.h"
#include "mscal_struc.h"
#include "l12_proto.h"

static int32_t nbands;

int rdstrideSDS(int32 fileID, char *sdsname, int32 stride[2], VOIDP array_data);
int read_glbl_attr(int32 sd_id, char *name, VOIDP ptr);
int crosscal_create(char *xcalfile, int32 *sd_id, mscalstr calstr, int32_t npixs);
int crosscal_readblocks(char *xcalfile, int32 *sd_id, int32_t *spix, int32_t *totalpixs, mscalstr *calstr);
int crosscal_writeblocks(char *xcalfile, int32 sd_id, int32_t *spix, mscalstr calstr, int32_t nfiles);


int crosscal_create(char *xcalfile, int32 *sd_id, mscalstr calstr, int32_t npixs)
{

    char title[255];
    l2prodstr p, *ptr;
    int j;
    

    sprintf(title, "%s Level-2 cross-calibration pixels", sensorName[calstr.sensorID]);

    
    *sd_id = SDstart(xcalfile, DFACC_CREATE);
    if (*sd_id == FAIL) {
       HY1C_out("-E- %s line %d: Could not create HDF file, %s .\n",__FILE__,__LINE__,xcalfile);
       return(HDF_FUNCTION_ERROR);
    }   
    
    /*                                                                  */
    /* Create the pixel SDSes                                           */
    /* ---------------------------------------------------------------- */
    /*                                                                  */
    idDS ds_id;
    ds_id.fid = *sd_id;
    ds_id.fftype = FMT_L3BIN;
    PTB( SetChrGA(ds_id, "Title", title) );
    PTB( SetI32GA(ds_id, "sensorID\0",(int32)calstr.sensorID) );
    PTB( SetChrGA(ds_id, "Input Parameters", calstr.input_parms));
    
    PTB( CreateSDS(
    *sd_id,                                          /* file id         */
    "fileID",                                        /* short name      */
    "ID of the File",                                /* long name       */
    NULL,                                            /* standard name   */
    "dimensionless",                                 /* units           */
    0, 0,                                            /* valid range     */
    0,0,                                             /* slope           */
    DFNT_INT16,                                      /* HDF number type */
    1,                                               /* rank            */
    npixs, 1, 1,                                     /* dimension sizes */
    "Number of Pixels", NULL, NULL                   /* dimension names */
    ) );

    PTB( CreateSDS(
    *sd_id,                                          /* file id         */
    "year",                                          /* short name      */
    "Pixel year",                                    /* long name       */
    NULL,                                            /* standard name   */
    "years",                                         /* units           */
    1996, 2038,                                      /* valid range     */
    0,0,                                             /* slope           */
    DFNT_INT16,                                      /* HDF number type */
    1,                                               /* rank            */
    npixs, 1, 1,                                     /* dimension sizes */
    "Number of Pixels", NULL, NULL                   /* dimension names */
    ) );

    PTB( CreateSDS(
    *sd_id,                                          /* file id         */
    "day",                                           /* short name      */
    "Pixel day of year",                             /* long name       */
    NULL,                                            /* standard name   */
     "days",                                         /* units           */
    0,366,                                           /* valid range     */
    0,0,                                             /* slope           */
    DFNT_INT16,                                      /* HDF number type */
    1,                                               /* rank            */
    npixs, 1, 1,                                     /* dimension sizes */
    "Number of Pixels", NULL, NULL                   /* dimension names */
    ) );

    PTB( CreateSDS(
    *sd_id,                                          /* file id         */
    "msec",                                          /* short name      */
    "Pixel time, milliseconds of day",               /* long name       */
    NULL,                                            /* standard name   */
    "milliseconds",                                  /* units           */
    0,0,                                             /* valid range     */
    0,0,                                             /* slope           */
    DFNT_INT32,                                      /* HDF number type */
    1,                                               /* rank            */
    npixs, 1, 1,                                     /* dimension sizes */
    "Number of Pixels", NULL, NULL                   /* dimension names */
    ) );

    PTB( CreateSDS(
    *sd_id,                                          /* file id         */
    "iscan",                                         /* short name      */
    "Scan-line number for the pixel",                /* long name       */
    NULL,                                            /* standard name   */
    "numbers",                                       /* units           */
    0,0,                                             /* valid range     */
    0,0,                                             /* slope           */
    DFNT_INT16,                                      /* HDF number type */
    1,                                               /* rank            */
    npixs, 1, 1,                                     /* dimension sizes */
    "Number of Pixels", NULL, NULL                   /* dimension names */
    ) );

    PTB( CreateSDS(
    *sd_id,                                          /* file id         */
    "mside",                                         /* short name      */
    "Mirror side for the pixel",                     /* long name       */
    NULL,                                            /* standard name   */
    "digit",                                         /* units           */
    0,0,                                             /* valid range     */
    0,0,                                             /* slope           */
    DFNT_UINT8,                                      /* HDF number type */
    1,                                               /* rank            */
    npixs, 1, 1,                                     /* dimension sizes */
    "Number of Pixels", NULL, NULL                   /* dimension names */
    ) );

    PTB( CreateSDS(
    *sd_id,                                          /* file id         */
    "detnum",                                        /* short name      */
    "Detector number for the pixel",                 /* long name       */
    NULL,                                            /* standard name   */
    "number",                                        /* units           */
    0,0,                                             /* valid range     */
    0,0,                                             /* slope           */
    DFNT_UINT8,                                      /* HDF number type */
    1,                                               /* rank            */
    npixs, 1, 1,                                     /* dimension sizes */
    "Number of Pixels", NULL, NULL                   /* dimension names */
    ) );

    PTB( CreateSDS(
    *sd_id,                                          /* file id         */
    "pixnum",                                        /* short name      */
    "Pixel number within the scan-line",             /* long name       */
    NULL,                                            /* standard name   */
    "number",                                        /* units           */
    0,0,                                             /* valid range     */
    0,0,                                             /* slope           */
    DFNT_INT16,                                      /* HDF number type */
    1,                                               /* rank            */
    npixs, 1, 1,                                     /* dimension sizes */
    "Number of Pixels", NULL, NULL                   /* dimension names */
    ) );

    PTB( CreateSDS(
    *sd_id,                                          /* file id         */
    "lon",                                           /* short name      */
    "Pixel longitude",                               /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -180.0, 180.0,                                   /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    1,                                               /* rank            */
    npixs, 1, 1,                                     /* dimension sizes */
    "Number of Pixels", NULL, NULL                   /* dimension names */
    ) );

    PTB( CreateSDS(
    *sd_id,                                          /* file id         */
    "lat",                                           /* short name      */
    "Pixel latitude",                                /* long name       */
    NULL,                                            /* standard name   */
    "degrees",                                       /* units           */
    -90.0, 90.0,                                     /* valid range     */
    0,0,                                             /* slope           */
    DFNT_FLOAT32,                                    /* HDF number type */
    1,                                               /* rank            */
    npixs, 1, 1,                                     /* dimension sizes */
    "Number of Pixels", NULL, NULL                   /* dimension names */
    ) );

    
    /*                                                                  */
    /* Create the geophysical SDSes for the L2 xcalibration products    */
    /* ---------------------------------------------------------------- */
    /*                                                                  */
    ptr = &p;
    for (j=0; j<calstr.nprods; j++) {
    
	if (strcmp(calstr.l2prods[j], "l2_flags") == 0) ;
	else
    	if (strcmp(calstr.l2prods[j], "mside") == 0) ;
	else
    	if (strcmp(calstr.l2prods[j], "detnum") == 0) ;
	else
    	if (strcmp(calstr.l2prods[j], "pixnum") == 0) ;
	else {
	
            if ((ptr = get_l2prod_index((char *)calstr.l2prods[j],(int32)calstr.sensorID,
                   (int32)calstr.nbands,0,1,calstr.Lambda)) == NULL) {

                fHY1C_out(stderr,
                "-E- %s line %d: product index failure.\n",
                __FILE__,__LINE__);
                return(FATAL_ERROR);
            };

            PTB( CreateSDS(
            *sd_id,                                      /* file id         */
            calstr.l2prods[j],                           /* short name      */
            ptr->title,                                  /* long name       */
            ptr->standard_name,                          /* standard name   */
            ptr->units,                                  /* units           */
            ptr->min, ptr->max,                          /* valid range     */
            0.0, 0.0,                                    /* slope, offset   */
            DFNT_FLOAT32,                                /* HDF number type */
            1,                                           /* number of dims  */
            npixs, 1, 1,                                 /* dimension sizes */
            "Number of Pixels", NULL, NULL               /* dimension names */
            ) );
	}
    }


    return(LIFE_IS_GOOD);
    
}





int crosscal_append(char *xcalfile, mscalstr calstr)
{

    int32 sd_id_old, sd_id;
    char  command[2*FILENAME_MAX], tempname[FILENAME_MAX];
    int   status, exists;
    int32_t  spix_old, spix, total_old, i;
    mscalstr calstr_old;
   
    nbands = NBANDS;   /* calstr.nbands; */

    spix_old = -1;
    total_old = 0;
    calstr_old.data = NULL;
    if (crosscal_readblocks(xcalfile, &sd_id_old, &spix_old, &total_old, &calstr_old) != LIFE_IS_GOOD) {
    	HY1C_out("-E- %s line %d: Cannot read already existing SDS data %s\n",__FILE__,__LINE__,xcalfile);
    	return(HDF_FUNCTION_ERROR);
    }

    if (calstr_old.npixs == 0L) exists = 0; else {
    	exists = 1;
    	if (calstr_old.sensorID != calstr.sensorID) {
    	    HY1C_out("-E- %s line %d: Attempt to write into a file whose data come from a different sensor %s\n",__FILE__,__LINE__,sensorName[calstr_old.sensorID]);
	    free_calstr(calstr_old, 1);
	    status = SDend(sd_id_old);
            return(HDF_FUNCTION_ERROR);
    	}
    }
    	    
    if (exists) { 
   
/*
        p = (char *)name;
	strcpy(tempname, xcalfile);
	p = dirname(tempname);
	strcpy(tempname, p);
	strcat(tempname, "/temp.hdf");
*/
	strcpy(tempname, xcalfile);
        strcat(tempname, "_tmp.hdf");
	
	if (crosscal_create(tempname, &sd_id, calstr, total_old+calstr.npixs) != LIFE_IS_GOOD) {
    	    HY1C_out("-E- %s line %d: Cannot read already existing SDS data %s\n",__FILE__,__LINE__,tempname);
	    free_calstr(calstr_old, 1);
	    status = SDend(sd_id);
	    status = SDend(sd_id_old);
            return(HDF_FUNCTION_ERROR);
    	}
	
	spix = 0L;
		
	do {
	    
	    if (crosscal_writeblocks(tempname, sd_id, &spix, calstr_old, 0L) != LIFE_IS_GOOD) {
    	    	HY1C_out("-E- %s line %d: Cannot write existing SDS data into the temp file %s\n",__FILE__,__LINE__,tempname);
	    	free_calstr(calstr_old, 1);
	    	status = SDend(sd_id);
	        status = SDend(sd_id_old);
    	    	return(HDF_FUNCTION_ERROR);
	    }
	    
	    if (crosscal_readblocks(xcalfile, &sd_id_old, &spix_old, &total_old, &calstr_old) != LIFE_IS_GOOD) {
    	    	HY1C_out("-E- %s line %d: Cannot read already existing SDS data %s\n",__FILE__,__LINE__,xcalfile);
	    	free_calstr(calstr_old, 1);
	    	status = SDend(sd_id);
    	    	return(HDF_FUNCTION_ERROR);
	    }

    	} while (calstr_old.npixs > 0L);


	for (i=0; i<calstr.npixs; i++) calstr.fileID[i] += (int16)calstr_old.nfiles;
	
	if (crosscal_writeblocks(tempname, sd_id, &spix, calstr, calstr_old.nfiles) != LIFE_IS_GOOD) {
    	    HY1C_out("-E- %s line %d: Cannot write new SDS data into the temp file %s\n",__FILE__,__LINE__,tempname);
	    status = SDend(sd_id);
    	    return(HDF_FUNCTION_ERROR);
	}
	    
    	if (SDend(sd_id) == FAIL) {
            HY1C_out("-E- %s line %d: Could not close HDF file, %s.\n",__FILE__,__LINE__,tempname);
            return(HDF_FUNCTION_ERROR);
    	}
	
	sprintf(command, "mv %s %s", tempname, xcalfile);
	system(command);

    } else {
    
    	if (crosscal_create(xcalfile, &sd_id, calstr, calstr.npixs) != LIFE_IS_GOOD) {
     	    HY1C_out("-E- %s line %d: Could not open HDF file, %s .\n",__FILE__,__LINE__,xcalfile);
	    status = SDend(sd_id);
            return(HDF_FUNCTION_ERROR);
    	}
	
	spix = 0L;
	
	if (crosscal_writeblocks(xcalfile, sd_id, &spix, calstr, 0L) != LIFE_IS_GOOD) {
    	    HY1C_out("-E- %s line %d: Cannot write new SDS data into the file %s\n",__FILE__,__LINE__,xcalfile);
	    status = SDend(sd_id);
    	    return(HDF_FUNCTION_ERROR);
	}
	    
    	if (SDend(sd_id) == FAIL) {
            HY1C_out("-E- %s line %d: Could not close HDF file, %s.\n",__FILE__,__LINE__,xcalfile);
            return(HDF_FUNCTION_ERROR);
    	}
    }
    
    
    return(LIFE_IS_GOOD);
    
}    



/* ------------------------------------------------------         */
/* crosscal_read() - reads cross-calibration data and 	          */
/*          subsamples if needed       			          */
/*          calptr comprises the actual data for npixs pixels     */
/*          memory is allocated depending on the npixs in the file*/
/* ------------------------------------------------------         */

int crosscal_read(char *xcalfile, int32_t subsmpl, mscalstr *calstr)
{

    int32 sd_id, sds_id;
    int32 dim_sizes[H4_MAX_VAR_DIMS];
    int32 stride[2], n_datasets, n_file_attr, rank, num_type, attributes;
    intn  status;
    char  name[H4_MAX_NC_NAME];
    int32_t  nfiles=0, totalpix=0, i, length, l;
    char  input_parms[16384];
    
    
    if (subsmpl <= 1) subsmpl = 1;
    
    sd_id = SDstart(xcalfile, DFACC_RDONLY);
    if (sd_id == FAIL) {
       	HY1C_out("-E- %s line %d: Could not open HDF file, %s .\n",__FILE__,__LINE__,xcalfile);
       	return(HDF_FUNCTION_ERROR);
    }
    
    strcpy(name, "sensorID\0");
    read_glbl_attr(sd_id, name, (VOIDP) &(calstr->sensorID));
    strcpy(name, "Input Parameters\0");
    read_glbl_attr(sd_id, name, (VOIDP) input_parms);
    l = strlen(input_parms);
    if ((calstr->input_parms = (char *)malloc((l+1)*sizeof(char))) == NULL) {
        HY1C_out("-E- %s line %d: Error allocating memory to MScalmerge input parameter text.\n",__FILE__,__LINE__);
    	return(HDF_FUNCTION_ERROR);
    }
    strncpy(calstr->input_parms, input_parms, l);
    calstr->input_parms[l] = '\x0';
    
    nfiles = 0L;
    do {
        sprintf(name, "filename%d", nfiles);
	    
  	if (SDfindattr(sd_id,name) == FAIL) break; else ++nfiles;
	
    } while (1);

    
    if (nfiles == 0L) {
       	HY1C_out("-E- %s line %d: Could not find filename attributes, %s .\n",__FILE__,__LINE__,name);
       	return(HDF_FUNCTION_ERROR);
    }
    
    
    strcpy(name, "fileID\0");
    if (getDims(sd_id, name, dim_sizes) != 0) {
       	HY1C_out("-E- %s line %d: Could not read HDF file dimensions, %s .\n",__FILE__,__LINE__,xcalfile);
        status = SDend(sd_id);
       	return(HDF_FUNCTION_ERROR);
    }
    if (dim_sizes[0] < 0){
    	HY1C_out("-E- %s line %d: The rank of the requested parameter (%s) is incorrect\n",__FILE__,__LINE__,name);
        status = SDend(sd_id);
    	return(HDF_FUNCTION_ERROR);
    }
    if (dim_sizes[0] == SD_UNLIMITED) HY1C_out("  Dimension SD_UNLIMITED"); else totalpix = (int32_t)dim_sizes[0];
    
        

    for (i=0; i<nfiles; i++) {
        sprintf(name, "filename%d", i);
	    
	if (read_glbl_attr(sd_id, name, (VOIDP) calstr->filenames[i]) != SUCCESS) {
	    HY1C_out("-E- %s line %d: Error reading filename attributes (%s)\n",__FILE__,__LINE__,name);
            free_calstr(*calstr, 1);
	    status = SDend(sd_id);
    	    return(HDF_FUNCTION_ERROR);
	}
    }
    

    calstr->nprods = 1000;
    if ((calstr->l2prods = (prname *)malloc(calstr->nprods*sizeof(prname))) == NULL) {
        HY1C_out("-E- %s line %d: Error allocating memory to l2 product names.\n",__FILE__,__LINE__);
    	return(HDF_FUNCTION_ERROR);
    }
    calstr->nprods = 0;
    
    SDfileinfo(sd_id, &n_datasets, &n_file_attr);
    
    
    for (i=0; i<n_datasets; i++) {
    
    	sds_id = SDselect(sd_id, i);
	status = SDgetinfo(sds_id, name, &rank, dim_sizes, &num_type, &attributes);
	
	if (strcmp(name, "fileID") == 0) ;
	else 
	if (strcmp(name, "l2_flags") == 0) ;
	else 
	if (strcmp(name, "year") == 0) ;
	else 
	if (strcmp(name, "day") == 0) ;
	else 
	if (strcmp(name, "msec") == 0) ;
	else 
	if (strcmp(name, "iscan") == 0) ;
	else 
	if (strcmp(name, "mside") == 0) ;
	else 
	if (strcmp(name, "detnum") == 0) ;
	else 
	if (strcmp(name, "pixnum") == 0) ;
	else 
	if (strcmp(name, "lon") == 0) ;
	else 
	if (strcmp(name, "lat") == 0) ;
	else {
	    
	    strcpy(calstr->l2prods[calstr->nprods], name);
	    calstr->nprods++;
    	}
    }
    
    if ((calstr->l2prods = (prname *)realloc((void *)calstr->l2prods, calstr->nprods*sizeof(prname))) == NULL) {
        HY1C_out("-E- %s line %d: Error reallocating memory to l2 product names.\n",__FILE__,__LINE__);
    	return(HDF_FUNCTION_ERROR);
    }
    calstr->Lambda = NULL;

	        
    totalpix = (totalpix+subsmpl-1)/subsmpl;
    
    length = alloc_calstr(nfiles, totalpix, calstr);
	    
    if (subsmpl <= 1) {
    	PTB( rdSDS(sd_id,"fileID\0",0,0,0,0,(VOIDP)calstr->fileID) );
    	PTB( rdSDS(sd_id,"year\0",0,0,0,0,(VOIDP)calstr->year) );
    	PTB( rdSDS(sd_id,"day\0",0,0,0,0,(VOIDP)calstr->day) );
    	PTB( rdSDS(sd_id,"msec\0",0,0,0,0,(VOIDP)calstr->msec) );
    	PTB( rdSDS(sd_id,"iscan\0",0,0,0,0,(VOIDP)calstr->iscan) );
    	PTB( rdSDS(sd_id,"mside\0",0,0,0,0,(VOIDP)calstr->mside) );
    	PTB( rdSDS(sd_id,"detnum\0",0,0,0,0,(VOIDP)calstr->detnum) );
    	PTB( rdSDS(sd_id,"pixnum\0",0,0,0,0,(VOIDP)calstr->pixnum) );
    	PTB( rdSDS(sd_id,"lon\0",0,0,0,0,(VOIDP)calstr->lon) );
    	PTB( rdSDS(sd_id,"lat\0",0,0,0,0,(VOIDP)calstr->lat) );
	   
	l = 0;
	for (i=0; i<calstr->nprods; i++) {
    	    if (strcmp(calstr->l2prods[i], "l2_flags") == 0) ;
    	    else
    	    if (strcmp(calstr->l2prods[i], "mside") == 0) ;
    	    else
    	    if (strcmp(calstr->l2prods[i], "detnum") == 0) ;
    	    else
    	    if (strcmp(calstr->l2prods[i], "pixnum") == 0) ;
    	    else {
	        PTB( rdSDS(sd_id,calstr->l2prods[i],0,0,0,0,(VOIDP)&(calstr->ddata[l*totalpix])) );
		++l;
	    }
	}
    } else {
    
	stride[0] = subsmpl;
	stride[1] = 1;
	PTB( rdstrideSDS(sd_id,"fileID\0",stride,(VOIDP)calstr->fileID) );
	PTB( rdstrideSDS(sd_id,"year\0",stride,(VOIDP)calstr->year) );
    	PTB( rdstrideSDS(sd_id,"day\0",stride,(VOIDP)calstr->day) );
    	PTB( rdstrideSDS(sd_id,"msec\0",stride,(VOIDP)calstr->msec) );
    	PTB( rdstrideSDS(sd_id,"iscan\0",stride,(VOIDP)calstr->iscan) );
    	PTB( rdstrideSDS(sd_id,"mside\0",stride,(VOIDP)calstr->mside) );
    	PTB( rdstrideSDS(sd_id,"detnum\0",stride,(VOIDP)calstr->detnum) );
    	PTB( rdstrideSDS(sd_id,"pixnum\0",stride,(VOIDP)calstr->pixnum) );
    	PTB( rdstrideSDS(sd_id,"lon\0",stride,(VOIDP)calstr->lon) );
    	PTB( rdstrideSDS(sd_id,"lat\0",stride,(VOIDP)calstr->lat) );
	   
	l = 0;
	for (i=0; i<calstr->nprods; i++) {
    	    if (strcmp(calstr->l2prods[i], "l2_flags") == 0) ;
    	    else
    	    if (strcmp(calstr->l2prods[i], "mside") == 0) ;
    	    else
    	    if (strcmp(calstr->l2prods[i], "detnum") == 0) ;
    	    else
    	    if (strcmp(calstr->l2prods[i], "pixnum") == 0) ;
    	    else {
	        PTB( rdstrideSDS(sd_id,calstr->l2prods[i],stride,(VOIDP)&(calstr->ddata[l*totalpix])) );
		++l;
	    }
	}
    }
    
    if (SDend(sd_id) == FAIL) {
       HY1C_out("-E- %s line %d: Could not close HDF file, %s .\n",__FILE__,__LINE__,xcalfile);
       return(HDF_FUNCTION_ERROR);
    }
    
    return(LIFE_IS_GOOD);
    
}    




/* ------------------------------------------------------         */
/* crosscal_npixs() - checks number of pixels contained in the    */
/*	    cross-calibration hdf file          	          */
/* ------------------------------------------------------         */

int crosscal_npixs(char *xcalfile, int32_t subsmpl, int32_t *sensorID, int32_t *npixs, int32_t *ngranuls)
{

    int32 sd_id;
    int32 dim_sizes[H4_MAX_VAR_DIMS];
    intn  status;
    char  name[H4_MAX_NC_NAME];
    int32_t  totalpix=0, totalfiles=0;
    
    
    if (subsmpl <= 1) subsmpl = 1;    
    
    sd_id = SDstart(xcalfile, DFACC_RDONLY);
    if (sd_id == FAIL) {
       	HY1C_out("-E- %s line %d: Could not open HDF file, %s .\n",__FILE__,__LINE__,xcalfile);
       	return(HDF_FUNCTION_ERROR);
    }
    
    strcpy(name, "sensorID\0");
    read_glbl_attr(sd_id, name, (VOIDP)sensorID);    
    
    totalfiles = 0L;
    do {
        sprintf(name, "filename%d", totalfiles);
	    
  	if (SDfindattr(sd_id, name) == FAIL) break; else ++totalfiles;
	
    } while (1);

    strcpy(name, "fileID\0");
    if (getDims(sd_id, name, dim_sizes) != 0) {
       	HY1C_out("-E- %s line %d: Could not read HDF file dimensions, %s .\n",__FILE__,__LINE__,xcalfile);
        status = SDend(sd_id);
       	return(HDF_FUNCTION_ERROR);
    }
    if (dim_sizes[0] < 0){
    	HY1C_out("-E- %s line %d: The rank of the requested parameter (%s) is incorrect\n",__FILE__,__LINE__,name);
        status = SDend(sd_id);
    	return(-1);
    }
    totalpix = (int32_t)dim_sizes[0];
    
    if (subsmpl > 1) totalpix = (totalpix+subsmpl-1)/subsmpl;
       
    if (SDend(sd_id) == FAIL) {
       HY1C_out("-E- %s line %d: Could not close HDF file, %s .\n",__FILE__,__LINE__,xcalfile);
       return(HDF_FUNCTION_ERROR);
    }
    
    *npixs = totalpix;
    *ngranuls = totalfiles;
        
    return(LIFE_IS_GOOD);
    
}    



/* ------------------------------------------------------         */
/* crosscal_readblocks() - reads consecutive cross-calibration    */ 	          
/*          blocks of data of the size of 400,000 pixels       	  */		        
/*          calptr comprises the actual data for the 400,0000     */
/*          pixels, memory is allocated with the first use        */
/* ------------------------------------------------------         */

int crosscal_readblocks(char *xcalfile, int32 *sd_id, int32_t *spix, int32_t *totalpixs, mscalstr *calstr)
{

    int32 dim_sizes[H4_MAX_VAR_DIMS], sds_id;
    intn  status;
    char  name[H4_MAX_NC_NAME];
    int32_t  nfiles=0, i, length, l;
    int32 n_datasets, n_file_attr, rank, num_type, attributes;
    char  input_parms[16384];
    
    
    if (*spix >= *totalpixs) {
    	if (SDend(*sd_id) == FAIL) {
            HY1C_out("-E- %s line %d: Could not close HDF file, %s .\n",__FILE__,__LINE__,xcalfile);
            return(HDF_FUNCTION_ERROR);
    	}
    	free_calstr(*calstr, 1);
        calstr->npixs = 0L;
	return(LIFE_IS_GOOD);
    }
	
    
    if (*spix < 0L) {
    
        calstr->npixs = 0L;
	calstr->nfiles = 0L;
	
    	*sd_id = SDstart(xcalfile, DFACC_RDONLY);
    	if (*sd_id == FAIL) {
	    return(LIFE_IS_GOOD);
    	}
    
    	strcpy(name, "sensorID\0");
    	if (read_glbl_attr(*sd_id, name, (VOIDP) &(calstr->sensorID)) != LIFE_IS_GOOD) {
       	    HY1C_out("-E- %s line %d: Could not read HDF sensor attribute, %s .\n",__FILE__,__LINE__,xcalfile);
            status = SDend(*sd_id);
       	    return(HDF_FUNCTION_ERROR);
    	}
	
    	strcpy(name, "Input Parameters\0");
    	read_glbl_attr(*sd_id, name, (VOIDP) input_parms);
    	length = strlen(input_parms);
    	if ((calstr->input_parms = (char *)malloc((length+1)*sizeof(char))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating memory to MScalmerge input parameter text.\n",__FILE__,__LINE__);
            status = SDend(*sd_id);
       	    return(HDF_FUNCTION_ERROR);
     	}
    	strncpy(calstr->input_parms, input_parms, length);
    	calstr->input_parms[length] = '\x0';
    
    
    	nfiles = 0L;
    	do {
            sprintf(name, "filename%d", nfiles);
	    
  	    if (SDfindattr(*sd_id,name) == FAIL) break; else ++nfiles;
	
    	} while (1);

    
    	strcpy(name, "fileID\0");
    	if (getDims(*sd_id, name, dim_sizes) != 0) {
       	    HY1C_out("-E- %s line %d: Could not read HDF file dimensions, %s .\n",__FILE__,__LINE__,xcalfile);
            status = SDend(*sd_id);
       	    return(HDF_FUNCTION_ERROR);
    	}
    	if (dim_sizes[0] < 0) {
    	    HY1C_out("-E- %s line %d: The rank of the requested parameter (%s) is incorrect\n",__FILE__,__LINE__,name);
            status = SDend(*sd_id);
    	    return(HDF_FUNCTION_ERROR);
     	}
    	*totalpixs = (int32_t)dim_sizes[0];
	
    	if (nfiles == 0L || *totalpixs == 0L) {
       	    HY1C_out("-E- %s line %d: There are no data in the file, %s .\n",__FILE__,__LINE__,name);
	    status = SDend(*sd_id);
      	    return(LIFE_IS_GOOD);
    	}

    	calstr->nprods = 1000;
    	if ((calstr->l2prods = (prname *)malloc(calstr->nprods*sizeof(prname))) == NULL) {
            HY1C_out("-E- %s line %d: Error allocating memory to l2 product names.\n",__FILE__,__LINE__);
            exit(FATAL_ERROR);
    	}
    	calstr->nprods = 0;
	
	
        SDfileinfo(*sd_id, &n_datasets, &n_file_attr);
    
    
    	for (i=0; i<n_datasets; i++) {
    
    	    sds_id = SDselect(*sd_id, i);
	    status = SDgetinfo(sds_id, name, &rank, dim_sizes, &num_type, &attributes);
	
	    if (strcmp(name, "fileID") == 0) ;
	    else 
	    if (strcmp(name, "l2_flags") == 0) ;
	    else 
	    if (strcmp(name, "year") == 0) ;
	    else 
	    if (strcmp(name, "day") == 0) ;
	    else 
	    if (strcmp(name, "msec") == 0) ;
	    else 
	    if (strcmp(name, "iscan") == 0) ;
	    else 
	    if (strcmp(name, "mside") == 0) ;
	    else 
	    if (strcmp(name, "detnum") == 0) ;
	    else 
	    if (strcmp(name, "pixnum") == 0) ;
	    else 
	    if (strcmp(name, "lon") == 0) ;
	    else 
	    if (strcmp(name, "lat") == 0) ;
	    else {
	    
	        strcpy(calstr->l2prods[calstr->nprods], name);
	        calstr->nprods++;
    	    }
        }
    
        if ((calstr->l2prods = (prname *)realloc((void *)calstr->l2prods, calstr->nprods*sizeof(prname))) == NULL) {
            HY1C_out("-E- %s line %d: Error reallocating memory to l2 product names.\n",__FILE__,__LINE__);
    	    return(HDF_FUNCTION_ERROR);
        }
        calstr->Lambda = NULL;

	
    	if (*totalpixs > 1000000) 
    	    length = alloc_calstr(nfiles, 1000000L, calstr);
	else
    	    length = alloc_calstr(nfiles, *totalpixs, calstr);
        

    	for (i=0; i<nfiles; i++) {
            sprintf(name, "filename%d", i);
	    
	    if (read_glbl_attr(*sd_id, name, (VOIDP) calstr->filenames[i]) != SUCCESS) {
	    	HY1C_out("-E- %s line %d: Error reading filename attributes (%s)\n",__FILE__,__LINE__,name);
            	status = SDend(*sd_id);
    		free_calstr(*calstr, 1);
        	calstr->npixs = 0L;
		calstr->nfiles = 0L;
    	    	return(HDF_FUNCTION_ERROR);
	    }
      	}
	
	*spix = 0L;
	
    } else {
    	
	if (*spix+calstr->npixs > *totalpixs) calstr->npixs = *totalpixs - *spix;
    }   
	    
    PTB( rdSDS(*sd_id,"fileID\0",(int32)*spix,0,calstr->npixs,0,(VOIDP)calstr->fileID) );
    PTB( rdSDS(*sd_id,"year\0",(int32)*spix,0,calstr->npixs,0,(VOIDP)calstr->year) );
    PTB( rdSDS(*sd_id,"day\0",(int32)*spix,0,calstr->npixs,0,(VOIDP)calstr->day) );
    PTB( rdSDS(*sd_id,"msec\0",(int32)*spix,0,calstr->npixs,0,(VOIDP)calstr->msec) );
    PTB( rdSDS(*sd_id,"iscan\0",(int32)*spix,0,calstr->npixs,0,(VOIDP)calstr->iscan) );
    PTB( rdSDS(*sd_id,"mside\0",(int32)*spix,0,calstr->npixs,0,(VOIDP)calstr->mside) );
    PTB( rdSDS(*sd_id,"detnum\0",(int32)*spix,0,calstr->npixs,0,(VOIDP)calstr->detnum) );
    PTB( rdSDS(*sd_id,"pixnum\0",(int32)*spix,0,calstr->npixs,0,(VOIDP)calstr->pixnum) );
    PTB( rdSDS(*sd_id,"lon\0",(int32)*spix,0,calstr->npixs,0,(VOIDP)calstr->lon) );
    PTB( rdSDS(*sd_id,"lat\0",(int32)*spix,0,calstr->npixs,0,(VOIDP)calstr->lat) );
	   
    l = 0;
    for (i=0; i<calstr->nprods; i++) {
    	if (strcmp(calstr->l2prods[i], "l2_flags") == 0) ;
    	else
    	if (strcmp(calstr->l2prods[i], "mside") == 0) ;
    	else
    	if (strcmp(calstr->l2prods[i], "detnum") == 0) ;
    	else
    	if (strcmp(calstr->l2prods[i], "pixnum") == 0) ;
    	else {
	    PTB( rdSDS(*sd_id,calstr->l2prods[i],(int32)*spix,0,calstr->npixs,0,(VOIDP)&(calstr->ddata[l*calstr->npixs])) );
	    ++l;
	}
    }   

    
    *spix += calstr->npixs;
        
    return(LIFE_IS_GOOD);
    
}    





/* ------------------------------------------------------         */
/* crosscal_writeblocks() - writes consecutive cross-calibration  */ 	          
/*          blocks of data                              	  */		        
/*          calptr comprises the actual data for the 400,0000     */
/*          pixels, memory is allocated with the first use        */
/* ------------------------------------------------------         */

int crosscal_writeblocks(char *xcalfile, int32 sd_id, int32_t *spix, mscalstr calstr, int32_t nfiles)
{

    char  name[H4_MAX_NC_NAME];
    int32_t  i, j, l;

    idDS ds_id;
    ds_id.fid = sd_id;
    ds_id.fftype = FMT_L3BIN;

    if (nfiles < 0) nfiles = 0L;	
    
    if (*spix <= 0 || nfiles > 0) {
    
    	for (i=0; i<calstr.nfiles; i++) {
    	
    	   j = i + nfiles;
           sprintf(name, "filename%d", j);
	    
    	   PTB( SetChrGA(ds_id, name, (char *)calstr.filenames[i]) );

	}
	    
    }

    PTB( sd_writedata(sd_id,"fileID\0",(VOIDP)calstr.fileID,(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
    PTB( sd_writedata(sd_id,"year\0",(VOIDP)calstr.year,(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
    PTB( sd_writedata(sd_id,"day\0",(VOIDP)calstr.day,(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
    PTB( sd_writedata(sd_id,"msec\0",(VOIDP)calstr.msec,(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
    PTB( sd_writedata(sd_id,"iscan\0",(VOIDP)calstr.iscan,(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
    PTB( sd_writedata(sd_id,"mside\0",(VOIDP)calstr.mside,(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
    PTB( sd_writedata(sd_id,"detnum\0",(VOIDP)calstr.detnum,(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
    PTB( sd_writedata(sd_id,"pixnum\0",(VOIDP)calstr.pixnum,(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
    PTB( sd_writedata(sd_id,"lon\0",(VOIDP)calstr.lon,(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
    PTB( sd_writedata(sd_id,"lat\0",(VOIDP)calstr.lat,(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
	   
    l = 0;
    for (i=0; i<calstr.nprods; i++) {
    	if (strcmp(calstr.l2prods[i], "l2_flags\0") == 0) ; 
	else
    	if (strcmp(calstr.l2prods[i], "mside\0") == 0) ; 
	else
    	if (strcmp(calstr.l2prods[i], "detnum\0") == 0) ;
	else
    	if (strcmp(calstr.l2prods[i], "pixnum\0") == 0) ;
	else {
	    PTB( sd_writedata(sd_id,calstr.l2prods[i],(VOIDP)(calstr.ddata+l*calstr.npixs),(int32)*spix,0,0,(int32)calstr.npixs,1,1) );
	    ++l;
	}
    }

    
    *spix += calstr.npixs;
    
    return(LIFE_IS_GOOD);
    
}    



/* ------------------------------------------------------  */
/* rdstrideSDS() - reads complete SDS (scientific data set)*/
/*           for a given stride and returns the data       */
/*                                                         */
/* ------------------------------------------------------  */
int rdstrideSDS(int32 fileID, char *sdsname, int32 stride[2], VOIDP array_data)
{
   int32 sds_id, status, numtype;
   int32 sds_index, rank, dims[H4_MAX_VAR_DIMS], nattrs;
   int32 start[2], edges[2];
   char  tmp_sdsname[H4_MAX_NC_NAME];


   /* Get the SDS index */
   sds_index = SDnametoindex(fileID,sdsname);

   if (sds_index < 0) {
     HY1C_out("-E- %s: SDS \"%s\" not found.\n", "rdSDS",sdsname);
     status = SDend(fileID);
     return(HDF_FUNCTION_ERROR);
   }


   /* Select the SDS */
   if ((sds_id = SDselect(fileID, sds_index)) == FAIL) {
     HY1C_out("-E- %s line %d: Could not read parameter %s in HDF file.\n",__FILE__,__LINE__,sdsname);
     status = SDend(fileID);
     return(HDF_FUNCTION_ERROR);
   }

   /* Get the rank and number type of the array */
   if ((SDgetinfo(sds_id, tmp_sdsname, &rank, dims, &numtype, &nattrs)) == FAIL) {
     HY1C_out("-E- %s line %d: Could not get info about parameter %s in HDF file.\n",__FILE__,__LINE__,sdsname);
     status = SDend(fileID);
     return(HDF_FUNCTION_ERROR);
   }


   /* Define the location, pattern and size of the data to read */
   /* set 1st dimension */
   start[0] = 0;
   start[1] = 0;
   edges[0] = dims[0];
   /* if rank > 1, set 2nd dimension */
   if (rank > 1) edges[1] = dims[1];

   /* Based on number type, call the corresponding wrapper
   for the HDF SDreaddata function */ 
   if ((SDreaddata(sds_id, start, stride, edges, array_data)) == FAIL) {
     HY1C_out("-E- %s line %d: Could not read data for parameter %s in HDF file.\n",__FILE__,__LINE__,sdsname);
     status = SDend(fileID);
     return(HDF_FUNCTION_ERROR);
   }
    
   /* Terminate access to the array */
   if ((SDendaccess(sds_id)) == FAIL) {
     HY1C_out("-E- %s line %d: Could not end access for parameter %s in HDF file.\n",__FILE__,__LINE__,sdsname);
     status = SDend(fileID);
     return(HDF_FUNCTION_ERROR);
   }

   return(LIFE_IS_GOOD);
   
}


int read_glbl_attr(int32 sd_id, char *name, VOIDP ptr) 
{                                           

  int32 index;
  
  if ((index = SDfindattr(sd_id,name)) == FAIL) return(HDF_FUNCTION_ERROR);

  if (SDreadattr(sd_id,index,(VOIDP)ptr)){ 
    HY1C_out("-E- %s line %d: Could not get global attribute, %s.\n", __FILE__,__LINE__,(name));
    return(HDF_FUNCTION_ERROR);
  }

  return(LIFE_IS_GOOD);
                                                         
}



int32_t alloc_calstr(int32_t nfiles, int32_t npixs, mscalstr *calstr)
{
    
    int32_t  len, l;
    unsigned char  *p;

    
    l = 5*sizeof(int16) + sizeof(int32) + 2*sizeof(uint8) + 2*sizeof(float);
    len = nfiles*sizeof(stname) + npixs*l;
    
    for (l=0; l<calstr->nprods; l++) {
    
    	if (strcmp(calstr->l2prods[l], "l2_flags") == 0) ;
	else
    	if (strcmp(calstr->l2prods[l], "mside") == 0) ;
        else
    	if (strcmp(calstr->l2prods[l], "detnum") == 0) ;
    	else
    	if (strcmp(calstr->l2prods[l], "pixnum") == 0) ;
    	else
	len += npixs*sizeof(float32);
    }
    
    if ((p = (unsigned char *) malloc(len)) == NULL) {
	HY1C_out("%s -Error: Cannot allocate memory to cross-calibration data\n",__FILE__);
   	exit(FATAL_ERROR);
    }
    calstr->nfiles = nfiles;
    calstr->npixs = npixs;
    calstr->data = p;
    calstr->filenames = (stname *) p; p += nfiles*sizeof(stname);
    calstr->fileID = 	(int16 	*) p; p += npixs*sizeof(int16);
    calstr->year = 	(int16 	*) p; p += npixs*sizeof(int16);
    calstr->day = 	(int16 	*) p; p += npixs*sizeof(int16);
    calstr->msec = 	(int32 	*) p; p += npixs*sizeof(int32);
    calstr->iscan = 	(int16 	*) p; p += npixs*sizeof(int16);
    calstr->mside = 	(uint8 	*) p; p += npixs*sizeof(uint8);
    calstr->detnum = 	(uint8 	*) p; p += npixs*sizeof(uint8);
    calstr->pixnum = 	(int16 	*) p; p += npixs*sizeof(int16);
    calstr->lon = 	(float 	*) p; p += npixs*sizeof(float);
    calstr->lat = 	(float 	*) p; p += npixs*sizeof(float);
    calstr->ddata = 	(float32 *) p;
    
    return(len);
    
}
   







void free_calstr(mscalstr calstr, int all)
{
     if (calstr.data != NULL) free(calstr.data);
     calstr.data = NULL;
     
     if (all) {
        if (calstr.Lambda != NULL) free(calstr.Lambda); calstr.Lambda = NULL;
	if (calstr.l2prods != NULL) free(calstr.l2prods); calstr.l2prods = NULL;
	if (calstr.input_parms != NULL) free(calstr.input_parms); calstr.input_parms = NULL;
     }
}


